# divide graphics window into 3 columns
par(mfrow=c(1,3))
# butted bar ends
par(lend=1)
 
# Poisson distribution when lambda = 1
plot(0:10, dpois(0:10, 1), type='h', lwd=6, xlab='Count', ylab='Probability')
mtext(side=3, line=.5, expression(lambda==1))
 
# Poisson distribution when lambda = 5
plot(0:15, dpois(0:15, 5), type='h', lwd=6, xlab='Count', ylab='Probability')
mtext(side=3, line=.5, expression(lambda==5))
 
# Poisson distribution when lambda = 15
plot(0:35, dpois(0:35, 15), type='h', lwd=6, xlab='Count', ylab='Probability')
mtext(side=3, line=.5, expression(lambda==15))
# return graphics window to normal state
par(mfrow=c(1,1))
 
# P(X = 16) for Poisson with lambda=15
dpois(16,15)
# normal approximation to Poisson probability
pnorm(16.5, 15, sqrt(15)) - pnorm(15.5, 15, sqrt(15))
 
# number of stems with various aphid counts
num.stems<-c(6,8,9,6,6,2,5,3,1,4)
# data frame of tabulated data
pois1 <- data.frame(aphids=0:9, counts=num.stems)
pois1
 
# construct raw data from tabulated data
aphid.data <- rep(pois1$aphids, pois1$counts)
aphid.data
length(aphid.data)
sum(num.stems)
 
# display the counts
barplot(num.stems)
# display the counts with labels
names(num.stems) <- 0:9
barplot(num.stems)
 
# Poisson probabilities of the data if lambda = 1
dpois(aphid.data, lambda=1)
# log probabilities
log(dpois(aphid.data, lambda=1))
# log-likelihood of the data when lambda = 1
sum(log(dpois(aphid.data, lambda=1)))
 
# log-likelihood of the data when lambda = 2
sum(log(dpois(aphid.data, lambda=2)))
 
# log-likelihood of the data when lambda = 3
sum(log(dpois(aphid.data, lambda=3)))
 
# function to calculate Poisson log-likelihood for aphids data
poisson.LL <- function(lambda) sum(log(dpois(aphid.data, lambda)))
 
# log-likelihood at lambda=3
poisson.LL(3)
 
# function gives incorrect answer with vector input
poisson.LL(seq(1,6,.1))
 
# investigating what went wrong
poisson.LL(1:2)
dpois(aphid.data,1)
dpois(aphid.data,2)
dpois(aphid.data,1:2)
 
# correct approach: use sapply
sapply(seq(1,6,.1), function(lambda) poisson.LL(lambda))
 
#plot the log-likelihood and look for maximum
plot(seq(1,6,.1), sapply(seq(1,6,.1), function(lambda) poisson.LL(lambda)), type='l', xlab=expression(lambda), ylab='log-likelihood')
# for Poisson MLE of lambda is the mean
mean(aphid.data)
 
### numerical optimization using R ####
?nlm
# nlm needs a function to minimize
pois.negloglike <- function(lambda) -poisson.LL(lambda)
# provide function and initial guess for estimate
out.pois <- nlm(pois.negloglike, 2)
# examine output
out.pois
# refit model to return Hessian
out.pois <- nlm(pois.negloglike, 2, hessian=T)
out.pois
# standard error of the mean
sqrt(1/out.pois$hessian)
 
#Wald confidence intervals
out.pois$estimate - 1.96*sqrt(1/out.pois$hessian)
out.pois$estimate + 1.96*sqrt(1/out.pois$hessian)